Coastal  and  Hydraulics  Laboratory  erdc/chltr-13-16 


US  Army  Corps 
of  Engineers® 

Engineer  Research  and 
Development  Center 


INNOVATIVE  SOLUTIONS 
for  a  safer,  better  world 


Military  Engineering 

Adaptive  Hydraulics/Hydrology  (AdH)  Pilot 
Point  Specification 

Guidelines  for  Solving  3D  Groundwater  Problems  Utilizing  Pilot  Points 

Kevin  D.  Winters  November  2013 


Approved  for  public  release;  distribution  is  unlimited. 


The  US  Army  Engineer  Research  and  Development  Center  (ERDC)  solves  the 
nation’s  toughest  engineering  and  environmental  challenges.  ERDC  develops  innovative 
solutions  in  civil  and  military  engineering,  geospatial  sciences,  water  resources,  and 
environmental  sciences  for  the  Army,  the  Department  of  Defense,  civilian  agencies,  and 
our  nation’s  public  good.  Find  out  more  at  vww.erdc.usace.armv.mil. 

To  search  for  other  technical  reports  published  by  ERDC,  visit  the  ERDC  online  library 
at  http://acwc.sdp.sirsi. net/client/default. 


Military  Engineering 


ERDC/CHL  TR-13-16 
November  2013 


Adaptive  Hydraulics/Hydrology  (AdH)  Pilot  Point 
Specification 

Guidelines  for  Solving  3D  Groundwater  Problems  Utilizing  Pilot  Points 


Kevin  D.  Winters 

Coastal  and  Hydraulics  Laboratory 

US  Army  Engineer  Research  and  Development  Center 

3909  Halls  Ferry  Road 

Vicksburg,  MS  39180-6199 


Final  report 

Approved  for  public  release;  distribution  is  unlimited. 


Prepared  for  US  Army  Corps  of  Engineers 

Washington,  DC  20314-1000 


ERDC/CHL  TR-13-16 


ii 


Abstract 

Guidelines  are  presented  herein  for  using  the  US  Army  Corps  of  Engineers 
(USACE)  Adaptive  Hydraulics/Hydrology  (AdH)  modeling  software  to 
model  three-dimensional  groundwater  problems  with  constituent  or  heat 
transport  utilizing  pilot  point  specification.  Pilot  point  specification  is  an 
auxiliary  module  intended  to  be  a  flexible  method  to  specify  spatially- 
varying  parameters  that  supersede  the  traditional  uniform  parameters  in 
the  model.  Examples  of  such  parameters  are  hydraulic  conductivity, 
porosity,  and  mesh  refinement.  Spatial  variation  can  be  used  to  develop 
high-fidelity  computer  models.  This  document  contains  descriptions  of  the 
pilot  point  input  cards  and  examples. 

The  pilot  point  specification  module  is  currently  integrated  into  the  AdH 
Groundwater  code  (kernel  version),  but  can  be  extended  to  additional  AdH 
physics  modules  as  necessary.  Input  is  currently  manually  generated  with 
result  viewable  utilizing  the  open-source  ParaView  visualization  software. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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1  Introduction 

Often,  a  groundwater  model  starts  with  the  process  of  simplifying  geologic 
data  (e.g.,  bore  logs)  into  hydrogeologic  units  to  create  a  conceptual  geologic 
model.  Each  unit  is  a  zone  of  similar  soil  characteristics  that  affect 
groundwater  flow.  This  zonal  classification  or  zonation  process  generalizes  a 
complex  system  into  regions  of  homogeneous  soil,  limiting  heterogeneity  to 
differences  between  hydrogeologic  units.  To  include  realistic  heterogeneity 
in  the  geologic  model,  many  zones  with  the  proper  spatial  distribution  may 
be  necessary.  Each  additional  zone  may  represent  the  same  basic  soil  type 
but  with  slightly  different  characteristics.  Once  a  conceptual  geologic  model 
is  finalized,  it  informs  the  creation  of  the  computational  domain  with  a 
partition  or  region  for  each  zone.  Each  domain  partition  is  assigned  a  set  of 
material  properties  endowing  the  domain  with  material  regions,  each  with 
uniform  characteristics. 

The  inclusion  of  spatially  varying  material  properties  with  smooth  transi¬ 
tions  is  not  practical  with  zonal  classification.  To  overcome  this  limitation, 
the  finite-element  software  framework  AdH  (Berger  and  Howington  2002; 
Pettway  et  al.  2010)  now  allows  specific  material  properties  to  vary 
smoothly  through  the  domain  with  the  use  of  pilot  points.  The  material 
region,  or  more  specifically,  the  partition  assigned  the  material  property, 
acts  as  the  ultimate  extent  of  the  variation.  In  this  way,  pilot  point  specifica¬ 
tion  complements  zonal  classification  by  permitting  heterogeneity  to  be 
defined  within  a  partition. 

Pilot  point  specification  consists  of  one  or  more  groups  of  the  following 
three  components: 

•  set  of  spatial  coordinates  with  values 

•  interpolation  method 

•  parameter  association. 

A  pilot  point  is  a  coordinate  location  with  a  given  value  or  values.  A 
collection  of  pilot  points  includes  enough  information  to  sufficiently 
describe  a  parameter  distribution  field  in  space  to  approximate  an 
unknown  physical  parameter.  Given  an  interpolation  method,  the 
provided  point  data  will  create  a  spatial  distribution  of  the  unknown 
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physical  parameter  on  the  computational  domain.  Once  associated  with  a 
material  property,  the  pilot  point  group  replaces  the  standard  uniform 
value  of  the  material  property;  a  value  is  now  specified  at  the  level  of  an 
individual  element  rather  than  the  material  region.  Each  pilot  point  group 
describes  a  unique  parameter  distribution  field,  such  as  the  variability  of 
hydraulic  conductivity  in  a  soil. 

While  pilot  point  specification  could  completely  replace  zonal 
classification,  by  declaring  the  entire  domain  as  a  single  partition  and 
describing  soil  properties  that  range  across  major  soil  types,  this  is,  in 
general,  not  advisable.  Not  all  material  properties  are  compatible  with 
pilot  point  specification  such  that  completely  different  soils  are  not 
possible.  For  example,  water  retention  curves  cannot  be  representative  of 
both  sand  and  clay  soils.  It  is  better  to  use  standard  zonal  classification  to 
describe  major  heterogeneity  (between  soil  types)  and  pilot  point 
specification  for  minor  heterogeneity  (internal  to  a  soil  type). 

The  original  concept  of  pilot  point  specification  focused  on  the  spatial 
variation  of  basic  soil  characteristics.  The  method’s  usefulness  has  been 
extended  to  include  most  material  properties,  mesh  adaption  control,  and 
initial  conditions. 

The  following  sections  provide  a  brief  example  of  pilot  points  incorporated 
in  a  simulation  (Section  2)  and  discussions  of  pilot  point  methods 
implemented  by  AdH,  input  control  cards,  output  data  (Section  3),  and 
AdH  model  execution  (Section  4). 
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2  Example  Simulation 

A  fictional  box  simulation  is  presented  here  to  showcase  aspects  of  pilot 
point  specification.  The  intent  of  this  example  is  to  contrast  the  results  of  a 
base  simulation  (standard  AdH  input)  and  simulations  utilizing  pilot 
points.  Familiarity  with  the  AdH  model  is  assumed  (documentation  and 
examples  available  at  http://adh.usace.armv.mil/).  A  detailed  discussion  of  pilot 
point  methods  and  specification  is  provided  in  Section  3. 

2.1  Base  Simulation 

The  domain  is  a  90-m  by  60-m  rectangular  prism  with  a  flat,  inclined 
surface  sloping  from  a  height  of  15  m  to  12  m  at  the  opposite  longitudinal 
edge  as  shown  in  Figure  1.  The  domain  includes  three  partitions,  repre¬ 
senting  a  silt  soil  (Material  1)  overlying  a  clay  soil  (Material  2),  overlying  a 
sand  soil  (Material  3).  A  flow  field  was  induced  by  specifying  head  values  of 
13.5  m  and  11  m  on  the  opposing  longitudinal  vertical  faces.  These  boundary 
conditions  are  hydrostatic  with  all  remaining  faces  assigned  no-flow 
boundary  conditions.  An  extraction  well  that  is  screened  in  the  bottom 
hydrogeologic  unit,  Material  3,  was  located  two-thirds  down  the  longitu¬ 
dinal  (x-axis)  centerline.  The  simulation  was  defined  using  the  standard 
AdH  input  cards  for  groundwater  problems,  including  zonal  hydraulic 
conductivities,  and  run  to  an  equilibrium  state  (Figure  1).  Material  2  acts  as 
an  aquitard  suppressing  the  effects  of  the  extraction  well  from  Material  1 
(the  effects  cannot  be  seen  in  Figure  1).  Material  3  is  the  preferential 
pathway  for  flow  since  the  sand  is  most  permeable. 

2.2  Inclusion  of  Pilot  Points 

Next,  pilot  point  specifications  were  included  to  depict  the  variation  of 
hydraulic  conductivities  in  each  hydrogeologic  unit.  For  convenience, 
hydraulic  conductivity  scaling  values  instead  of  hydraulic  conductivity 
parameter  values  are  given  at  each  domain  corner  and  the  well  for  each 
material  type  (Section  3  details  the  available  options  for  specifying  hydraulic 
conductivity).  The  resulting  hydraulic  conductivity  values  at  the  domain 
corners  are  listed  in  Table  1.  Cardinal  directions  are  used  to  denote  the 
domain  corners  with  the  positive  y-axis  direction  aligned  with  North.  It  is 
noted  that  the  pilot  points  located  at  the  well  have  scalar  values  of  1.0  that 
reproduce  the  zonal  hydraulic  conductivity  values  used  in  the  base 
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Figure  1.  Example  computational  domain  showing  zonal  hydraulic  conductivities  (Kh  )  and 
pilot  point  locations  as  yellow  cylinders  (A),  and  the  computed  total  head  solution  (B). 


Table  1.  Specified  Hydraulic  Conductivities  for  Example  Simulation 


Material 

Horizontal  Hydraulic  Conductivities,  Kh  (m/day) 

Zonal 

Pilot  Points 

NW 

NE 

SE 

SW 

Well 

1  (silt) 

0.300 

0.0660 

0.276 

1.01 

0.645 

0.300 

2  (clay) 

1.00e-5 

5.40e-6 

2.27e-5 

1.60e-6 

5.32e-5 

1.00e-5 

3  (sand) 

5.00 

0.81 

1.74 

37.7 

15.9 

5.00 

simulation.  The  five  point  locations  are  depicted  in  Figure  l.  Each  material 
type  has  a  pilot  point  group  associated  with  the  hydraulic  conductivity 
tensor.  The  scaling  values  are  interpolated  to  the  elements  assigned  the 
respective  material  type  and  alter  the  hydraulic  conductivity  tensor  used  in 
the  system  of  groundwater  equations.  The  solution  will  then  be  based  on 
values  specific  to  each  element  as  well  as  the  non- varying  material 
properties. 

The  box  simulation  was  run  with  each  available  interpolation  method 
(natural  neighbor,  inverse-distance  weighted,  and  ordinary  kriging; 
described  in  Section  3)  using  the  same  pilot  point  values.  The  resulting 
hydraulic  conductivity  fields  generated  by  pilot  point  specification  are 
displayed  in  Figures  2  through  4,  contrasted  by  the  original  base  simulation 
zonal  values.  Pilot  point  specification  created  asymmetrical,  macroscopic 
effective  conductivity  zones  across  each  material  layer  (hydraulic  conduc¬ 
tivity  is  still  isotropic  locally).  Figures  5  and  6  show  the  resulting  head 
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Figure  2.  Comparison  of  the  horizontal  hydraulic  conductivity  (Kh)  of  Material  1  using  a  variety  of 
methods:  (A)  zonal,  (B)  nearest-neighbor,  (C)  inverse-distance  weighted,  and  (D)  ordinary  kriging. 


Figure  3.  Comparison  of  the  horizontal  hydraulic  conductivity  (Kh)  of  Material  2  using  a  variety  of 
methods:  (A)  zonal,  (B)  nearest-neighbor,  (C)  inverse-distance  weighted,  and  (D)  ordinary  kriging. 
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Figure  4.  Comparison  of  the  horizontal  hydraulic  conductivity  (Kh)  of  Material  3  using  a  variety  of 
methods:  (A)  zonal,  (B)  nearest-neighbor,  (C)  inverse-distance  weighted,  and  (D)  ordinary  kriging. 


Figure  5.  Comparison  of  the  total  head  solution  based  on  a  variety  of  methods  to  describe 
hydraulic  conductivity:  (A)  zonal,  (B)  nearest-neighbor,  (C)  inverse-distance  weighted,  and 

(D)  ordinary  kriging. 
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Figure  6.  Comparison  of  the  total  head  solution,  clipped  through  Material  3  to  show  influence 
of  extraction  well,  based  on  a  variety  of  methods  to  describe  hydraulic  conductivity:  (A)  zonal, 
(B)  nearest-neighbor,  (C)  inverse-distance  weighted,  and  (D)  ordinary  kriging. 


solutions  of  the  full  domain  and  the  domain  clipped  at  3.5  m  to  highlight  the 
influence  of  the  extraction  well,  respectively.  Looking  specifically  at  the 
lowest  hydrogeologic  unit,  Material  3,  the  permutations  generated  zones  of 
higher  and  lower  permeability  generally  parallel  to  the  macroscopic  flow 
field  (longitudinal  axis).  The  connectivity  of  the  similar  zones  and  the 
smoothness  of  variations  are  dependent  on  the  interpolation  method 
employed.  These  hydraulic  conductivity  zones  bend  the  head  contours  from 
the  original  symmetrical  solution  producing  differing  well-capture  areas  in 
the  confined  aquifer. 

2.3  Review 

While  this  example  shows  the  successful  inclusion  of  a  spatially  varying 
material  property,  it  is  also  provided  as  a  cautionary  tale.  It  is  easy  to 
append  a  standard  AdH  simulation  with  pilot  point  specification  and 
generate  complex  soil  characteristics,  but  it  may  not  be  beneficial.  The 
domain  is  split  by  an  aquitard  such  that  the  hydraulic  conductivity  of  the 
layers  may  be  inconsequential  beyond  their  relative  magnitudes  depending 
on  the  purpose  of  the  model.  For  example,  if  the  intent  is  to  verily  the 
extents  of  the  capture  zone  of  the  extraction  well,  then  pilot  points  only  for 
Material  3  could  be  added  after  verifying  that  it  is  indeed  within  a  confined 
aquifer.  On  the  other  hand,  pilot  points  could  be  added  only  for  Material  1 
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when  determining  the  location  of  the  water  table.  Additionally,  the 
permeability  of  the  clay,  Material  2,  is  almost  small  enough  to  be 
considered  a  solid,  computationally.  Though  this  problem  size  is  small, 
each  pilot  point  parameter  taxes  the  model  resources  by  requiring 
additional  memory  allocation,  computations,  and  logical  operations.  These 
effects  are  scaled  as  the  problem  size  increases.  Finally,  it  is  important  to 
utilize  the  appropriate  interpolation  method  to  describe  significant  spatial 
variation  of  the  distribution  field. 

Ideally,  pilot  point  specification  for  material  properties  would  be  used  in 
conjunction  with  parameter  estimation  methods  to  calibrate  results  with 
observed  data.  New  parameter  values  can  be  iteratively  replaced  in  the 
pilot  point  specification  without  the  need  to  edit  the  original  AdH  input 
definition  of  a  soil.  Thus,  pilot  points  serve  to  regularize  an  otherwise 
difficult  parameter  estimation  problem.  Pilot  point  specification  should  be 
included  only  after  the  groundwater  system  is  well  understood,  the  regular 
zonal  classified  simulation  is  run,  and  the  results  indicate  specific  areas  of 
the  domain  that  need  to  be  addressed  further. 
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3  Methods  and  Files 

To  assist  the  user,  this  section  combines  two  essential  discussions  of  pilot 
point  specification,  the  implementation  of  methods  within  AdH  and  the 
input  and  output  files.  The  desire  is  to  ensure  clarity  when  seeking  informa¬ 
tion  of  a  given  pilot  point  card.  File  cards  are  bolded  in  the  document  body 
for  ease  of  use. 

3.1  Input 

Since  pilot  point  specifications  supersede  the  standard  input,  the  location 
of  the  pilot  point  cards  was  designed  to  be  flexible.  The  default  location  to 
provide  pilot  point  cards  is  a  new  ASCII  text  file  given  the  simulation  base 
name  with  the  designation  “pp”  as  the  extension.  For  example,  if  the  AdH 
simulation  is  named  “my_sim”  (i.e.,  my_sim.bc)  then  the  default  location 
is  the  file  my_sim.pp.  Alternatively,  the  pilot  point  cards  may  be  provided 
in  files  otherwise  named  if  referenced  within  the  AdH  super  file  (Asup) 
with  the  PP  card  (Example  B  and  Example  C  in  the  Appendix).  This 
flexibility  permits  the  specification  of  pilot  points  within  the  standard  AdH 
input  file  (Abe)  or  across  multiple  files.  See  the  appendix  for  examples  of 
input  file  combinations. 

Currently,  the  pilot  point  input  must  be  manually  generated;  a  complete 
graphical  user  interface  (GUI)  to  generate  the  input  will  eventually  be 
available  in  the  Computational  Model  Builder  (CMB,  developed  by  Kitware, 
Inc.,  for  ERDC)  suite’s  ModelBuilder  tool.  Specified  pilot  point  information 
supersedes  original  information  only  if  the  operational  parameter  card  OP 
PP  is  inserted  in  the  input  file;  otherwise,  provided  pilot  point  information 
is  ignored,  and  the  original  information  retains  its  precedence. 

The  pilot  point  specification  cards,  listed  in  Table  2  and  described  in  the 
following  subsections,  may  appear  in  any  order  with  the  cards  forming  sets 
linked  together  by  unique  pilot  point  group  IDs.  Pilot  point  specification 
utilizes  IDs  that  are  not  restricted  to  the  one-based,  sequential  limitations  of 
other  AdH  IDs;  IDs  may  be  positive  or  negative  and  are  limited  only  by  the 
operating  system’s  definition  of  an  integer.  It  is  best  practice  to  choose 
predetermined  ranges  of  IDs  to  represent  parameters  for  easy  inclusion  into 
the  simulation.  For  example,  single  or  double  digit  IDs  could  represent 
specific  initial  conditions  while  larger  numbers  could  refer  to  specific 
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materials  by  incorporating  the  material  ID  (25010  and  25020  could  allude 
to  material  25).  In  this  way,  different  combinations  of  pilot  point 
specification  can  be  included  without  renumbering  the  input. 


Table  2.  Control  Card  Categories 


Card 


Description 


Operation  Parameters 


(Section  3.1.1  and  Table  31 


OP  PP 


Enable  Pilot  Point  Specification 


Pilot  Point  Specification 


(Section  3.1.21 


Association  Parameters 


(Section  3.1.2.1  and  Table  41 


PP  HOT 


Initial  Condition  Parameter 


PP  MP 


Material  Parameter 


Group  Properties 


(Section  3.1.2.2  and  Table  5) 


PP  LIM 


Interpolation  Limits 


PP  PT2 


2D  Pilot  Points 


PP  PT3 


3D  Pilot  Points 


PP  RAD 


Search  Radius 


PPTYP 


Interpolation  Method 


Kriging  Interpolation  Properties 


(Section  3.1.2.3  and  Table  6) 


PP  KRG 


Kriging  Information 


PP  VGC 


Variogram  Contributions 


PP  VGI 


Variogram  IDs 


PP  VGS 


Variogram  Sill 


PP  VGW 


Variogram  Weights 


PP  VG2 


2D  Variogram  Information 


PP  VG3 


Miscellaneous 


PP  DBG 


3D  Variogram  Information 


(Section  3.1.2.4  and  Table  7) 


Debug  Information 


Table  entries  are  hyperlinked. 


Comments  are  permitted  in  the  input  files  if  demarcated  with  a  preceding 
#  or  /;  all  text  after  the  delimiter  on  the  file  line  is  ignored.  Blank  lines  are 
also  permitted.  The  AdH  pilot  point  file  input  routines  will  validate  the 
cards  and  provide  information  and  error  messages  to  assist  correct  card 
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specification.  AdH  will  exit  after  listing  any  errors  and  before  pilot  point 
interpolation  and  normal  simulation  operation  occurs. 

3.1.1  Operation  Parameters 

The  problem  type  and  operational  methods  of  AdH  are  controlled  by  the 
operational  parameter  cards,  which  are  denoted  by  OP  card  type.  To 
utilize  the  pilot  point  specification  in  a  groundwater  and/or  heat  transport 
model,  the  OP  PP  card  (Table  3)  must  be  included  with  the  normal 
operational  parameter  cards.  This  card  is  the  flag  for  AdH  to  perform  the 
auxiliary  logical  operations  to  support  pilot  points.  AdH  performs  the 
standard  model  operation,  specified  by  the  normal  operation  parameters 
(e.g.,  OP  GW)  if  pilot  point  cards  (PP  card  type)  are  given  and  OP  PP 
card  is  excluded. 


Table  3.  Operation  Parameter  Cards 


OPPP 

ENABLE  PILOT  POINT  SPECIFICATION 

The  OP  PP  card  enables  the  operation  of  pilot  point  specification  and  must  be  included 
only  in  the  standard  AdH  input  file  (*.bc). 

Field 

Type 

Value 

Description 

1 

string 

OP 

Card  type 

2 

string 

PP 

Parameter 

3.1.2  Pilot  Point  Specification 

Pilot  point  specification  cards  are  identified  by  the  designation  PP  and  are 
sorted  into  subcategories  described  in  the  following  subsections.  There  will 
be  a  set  of  cards  for  each  pilot  point  group  included  in  the  model. 

3. 1.2.1  Association  Parameters 

These  cards  specify  the  spatially  varying  simulation  parameter  that  a  pilot 
point  group  represents.  A  pilot  point  group  may  only  be  associated  with  a 
single  parameter:  an  initial  condition  or  a  material  property. 

Pilot  point  groups  may  describe  any  of  the  available  initial  condition  types 
for  groundwater  (pressure  head  or  total  head  and  concentration)  and  heat 
transport  problems  (temperature)  by  the  PP  HOT  card.  The  specification  is 
applied  to  the  entire  domain  by  interpolating  values  at  all  node  locations. 
The  standard  hot  file  is  still  required  by  AdH;  therefore,  pilot  point 
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specifications  will  supersede  or  supplement  the  model’s  regular  initial 
conditions. 

AdH  material  properties  are  divided  into  two  sets:  global  (e.g.,  gravity,  MP 
G)  and  material-specific  (e.g.,  hydraulic  conductivity,  MP  K).  Pilot  point 
groups  may  describe  a  subset  of  groundwater  and  heat  problem  material 
specific  properties  that  are  suitable  for  spatial  interpolation  with  the  PP 
MP  card,  including  the  following: 

•  maximum  refinement  level 

•  refinement  tolerance 

•  hydraulic  conductivity 

•  porosity 

•  specific  storage 

•  dispersivity 

•  tortuosity 

•  molecular  diffusion 

•  retardation  coefficient. 

The  pilot  point  group  specification  is  confined  to  the  given  material’s 
region(s)  in  the  domain;  the  new  property  is  interpolated  at  the  centroid  of 
each  element  assigned  the  material.  Examples  of  an  unsuitable  material 
property  are  the  water  retention  curves  (pressure-relative  conductivity  and 
pressure-saturation)  since  X-Y  series  are  necessary  to  describe  these 
relationships,  even  if  van  Genuchten  parameters  are  used.  The  majority  of 
the  permitted  material  properties  are  single  real-data-type  values 
representing  characteristics  that  in  reality  are  spatially  heterogeneous  and 
therefore  are  a  perfect  fit  for  pilot  point  specification;  the  following  two 
cases  are  not  as  direct. 

The  key  ability  of  AdH,  the  adaption  of  the  mesh,  is  controlled  by 
tolerances  and  level  flags.  Mesh  elements  are  split,  or  refined,  when  the 
calculated  error  indicators  at  the  elements  are  greater  than  a  given 
refinement  tolerance  (e.g.,  MP  FRT),  and  the  elements’  refinement  levels 
are  less  than  the  given  maximum  (MP  ML)  where  the  level  of  an  element 
is  the  number  of  times  an  element  has  been  refined.  A  maximum 
refinement  level  of  zero  eliminates  adaption.  Mesh  elements  are  merged, 
or  unrefined,  when  the  calculated  error  indicators  are  significantly  smaller 
than  the  refinement  tolerance.  Pilot  point  specification  may  be  utilized  to 
control  mesh  adaption  within  a  material  region  by  spatially  varying  a 
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specific  problem  tolerance  and/or  the  maximum  level.  The  maximum  level 
parameter  is  a  discrete  quantity;  hence,  the  pilot  point  scheme’s 
interpolated  value  (a  real  data  type)  is  rounded  up  or  down  to  the  nearest 
integral  value  prior  to  its  substitution. 

The  standard  hydraulic  conductivity  card,  MP  K,  requires  six  real-data-type 
values  to  describe  the  second-order  symmetric  tensor  (Kxy  =  Kyx,  Kxz  =  Kzx, 
and  Kyz  =  Kzy ).  Tensor  interpolation  is  not  supported  by  AdH’s  pilot  point 
specification,  so  three  alternatives  are  provided  to  specify  heterogeneity. 
Option  l:  the  originally  specified  hydraulic  conductivity  tensor  may  be 
scaled  by  a  spatially  varying  factor.  All  tensor  components  are  multiplied  by 
the  same  factor.  Option  2:  the  tensor  may  be  superseded  by  two  separate 
spatially  varying  horizontal  and  vertical  conductivities  ( Kh  and  Kv, 
respectively)  where  their  respective  off-diagonal  tensor  components  are 
defaulted  to  zero.  If  Kh  is  given,  then  Kxx  =  Kyy  =  Kh  and  Kxy  =  Kyx  =  Kxz  = 
Kzx  =  Kyz  =  Kzy  =  O,  while  Kxz  is  unchanged  as  shown  in  Figure  7.  If  Kv  is 
given,  then  Kzz  =  Kv  and  Kxz  =  Kzx  =  Kyz  =  Kzy  =  O  with  the  remaining 
components  unchanged.  Option  3:  the  tensor  may  be  superseded  at  the 
individual  component  level  (Kxx,  Kxy,  Kxz,  Kyy,  Kyz,  Kzz).  Multiple  pilot 
point  groups  will  be  necessary  to  completely  specify  a  material’s  hydraulic 
conductivity  tensor.  These  three  options  are  mutually  exclusive  at  the 
material  level;  for  a  given  material,  specifying  a  scaling  factor,  horizontal 
conductivity,  and  Kzz  component  are  prohibited,  but  material  1  may  specify 
a  scaling  factor;  material  2,  horizontal  conductivity;  and  material  3,  an 
individual  component.  It  is  not  required  that  both  Kh  and  Kv,  and,  similarly, 
all  the  individual  components,  are  specified.  The  original  hydraulic 
conductivity  may  be  partially  superseded  by  pilot  point  specification. 

The  association  parameter  cards  are  mutually  exclusive  with  one  required 
for  each  pilot  point  group  (Table  4). 

3. 1.2. 2  Group  Properties 

These  cards  define  the  basic  information  of  a  pilot  point  group  (a  set  of 
points  with  values  representative  of  a  given  parameter)  and  the  method  to 
interpolate  parameter  values  to  the  domain. 

Pilot  point  specification  may  utilize  a  set  of  either  2D  or  3D  point  locations, 
given  by  the  PP  PT2  or  PP  PT3  cards,  respectively,  though  the  available 
interpolation  methods  are  limited  for  the  lower  dimension.  In  either  case 
(2D  or  3D),  a  list  of  labels,  coordinates,  and  parameter  values  is  required. 
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Figure  7.  Hydraulic  conductivity,  K,  tensor:  (A)  original;  altered  by  option  1  (B) 
scale  factoring  (in  red,  bold);  by  option  2  (C)  horizontal  conductivity  and  (D) 
vertical  conductivity;  and  by  option3  individual  components  (E)  Kxx  (red),  Kyy 
(green),  Kzz  (blue),  (F)  Kxy  (red),  Kxz  (green),  and  Kyz  (blue). 
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Table  4.  Association  Parameter  Cards 


PP  HOT 

INITIAL  CONDITION  PARAMETER 

The  PP  HOT  card  specifies  the  initial  domain  condition  a  pilot  point  group  represents  for 
hot  starting  the  simulation.  The  last  field  is  conditional  on  the  penultimate  field. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

HOT 

Parameter 

3 

string 

Initial  condition: 

IPH  Pressure  head  (GW) 

ITH  Total  head  (GW) 

IC  Constituent  (GW) 

IT  Temperature  (Heat) 

If  initial  condition  is  constituent  (field  3  is  equal  to  IC]: 

4 

int 

>0 

Constituent  ID 

PP  MP 

MATERIAL  PARAMETER 

The  PP  MP  card  specifies  the  material  property  a  pilot  point  group  represents.  The  last 
field  is  conditional  on  a  preceding  field.  It  is  noted  that  material  IDs  within  AdH  are  no 
longer  restricted  to  a  consecutive  series  beginning  with  1  but  can  be  any  integer  value 
except  the  C  language  macro  constant  INT_MAX  value  (found  in  limits. h;  actual  value 
depends  on  operating  system  and  library  implemention). 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

MP 

Parameter 
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3 

string 

Material  property: 

DF 

Molecular  diffusion  (GW) 

DPL 

Longitudinal  dispersivity  (GW  &  Heat) 

DPT 

Transverse  dispersivity  (GW  &  Heat) 

FRT 

Flow  refinement  tolerance  (GW) 

HRT 

Heat  refinement  tolerance  (Heat) 

KS 

Hydraulic  conductivity  scaling 

factor  (GW) 

KH 

Horizontal  hyd.  cond.  (GW) 

KV 

Vertical  hyd.  cond.  (GW) 

KXX 

Hyd.  cond.  tensor  XX  component  (GW) 

KXY 

Hyd.  cond.  tensor  XY  comp.  (GW) 

KXZ 

Hyd.  cond.  tensor  XZ  comp.  (GW) 

KYY 

Hyd.  cond.  tensor  YY  comp.  (GW) 

KYZ 

Hyd.  cond.  tensor  YZ  comp.  (GW) 

KZZ 

Hyd.  cond.  tensor  ZZ  comp.  (GW) 

ML 

Max.  level  of  refinement  (GW  &  Heat) 

POR 

Porosity  (GW  and  Heat) 

RD 

Retardation  (GW) 

SS 

Specific  storage  (GW) 

TOR 

Tortuosity  (GW  &  Heat) 

TRT 

Transport  refinement  tolerance  (GW) 

4 

int 

*  INT_MAX 

Material  ID 

If  material  property  is  molecular  diffusion,  retardation,  or  transport  refinement  tolerance 

(field  3  is  equal  to  DF,  RD,  or  TRT,  respectively): 

5 

Int 

>0 

Constituent  ID 

The  point  labels  are  ignored  by  AdH  but  included  to  assist  in  the  classifica¬ 
tion  and  visualization  of  points  in  other  software.  It  is  best  practice  to  export 
point  information  from  third-party  applications  such  as  Microsoft  Excel  or 
the  Groundwater  Modeling  System  (GMS). 

The  PP  RAD  card  apportions  the  subset  of  pilot  points  involved  in  any 
given  interpolation  with  a  search  radius  and  point  count  limits  while  the  PP 
LIM  card  defaults  or  restricts  the  interpolated  value.  As  shown  in  Figure  8, 
the  search  radius  must  be  chosen  wisely  as  it  directly  influences  interpola¬ 
tion  at  the  locations  of  interest.  If  the  radius  is  too  small,  no  pilot  points  may 
be  found,  and,  therefore,  interpolation  cannot  occur.  Additionally,  if  the 
radius  encircles  fewer  points  than  the  given  minimum  requirement  (e.g., 
specified  minimum  of  3  points,  and  the  green  circle  in  the  Figure  8  is  given), 
then  interpolation  will  not  occur.  In  the  case  of  a  large  subset  of  pilot  points 
where  the  count  exceeds  the  given  maximum  limit,  the  points  are  prioritized 
by  distance  from  the  location  of  interest  with  the  farthest  points  removed 
from  the  subset  until  the  count  is  equal  to  the  maximum  count  limit.  AdH 
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does  not  include  a  method  to  sample  pilot  points  by  quadrant  nor 
guarantees  that  the  subset  completely  surrounds  the  location  of  interest. 
The  distance  between  a  pilot  point  and  the  location  of  interest  is  computed 
by  Equations  l  and  2,  based  on  the  given  interpolation  scheme 
dimensionality. 


Figure  8.  Example  pilot  point  search  radii  for  locations  of  interest  (blue  point)  that 
encircle  zero,  one,  and  four  points  (read,  green,  and  blue  dashed  circles  respectively). 


3D: 


d  = 


2D  horizontal: 


d  = 


(1) 


(2) 


where: 
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d.  =  distance,  L 

xloi>  yiov  zioi  -  coordinates  of  the  location  of  interest,  L 
xpt>  yPt.  zpt  -  coordinates  of  the  pilot  point,  L. 

When  search  criteria  are  not  satisfied  and  interpolation  fails,  the  user- 
defined  default  is  assigned  at  the  location  of  interest.  Two  use  cases  must 
be  acknowledged  for  this  default  value.  First,  specifying  a  blatantly  invalid 
value  can  test  a  pilot  point  group’s  interpolation  scheme  and  ensure 
complete  coverage  of  the  domain  (or  material  region).  The  value  may 
cause  premature  termination  of  AdH.  Second,  if  desired  spatial  variance  is 
limited,  the  default  value  can  provide  the  standard  parameter  value  to  be 
augmented  by  interpolation.  In  other  words,  the  default  value  can  be 
either  a  flag  for  gaps  or  an  intentional  fill  value. 

After  successful  interpolation,  the  computed  value  is  compared  with  given 
interpolation  limits  and  restricted  to  the  range  as  necessary.  For  example, 
given  the  desired  range  of  to  to  20,  an  interpolated  value  of  21  will  be 
revised  down  to  20. 

The  interpolation  method  is  declared  using  the  PP  TYP  card  and  consists 
of  two  components:  type  and  dimensionality.  AdH  currently  provides 
three  interpolation  types  that  do  not  depend  on  a  priori  relationships  of 
data  points  (e.g.,  Delaunay  triangulation  or  Voronoi  diagram): 

•  nearest-neighbor 

•  modified  Shepard’s  method  inverse-distance  weighted  (IDW) 

•  ordinary  kriging 

and  two  aforementioned  dimensionalities: 

•  3D 

•  2D  horizontal 

for  a  combinatorial  total  of  six  schemes.  The  dimensionality  component 
restricts  which  coordinates  are  used  during  computation,  specifically  in 
regard  to  distances.  Three-dimensionality  allows  all  coordinates  to  be 
included;  2D  horizontal  dimensionality  ignores  all  z-coordinate  (elevation) 
values  by  defaulting  them  to  zero  prior  to  computations.  The  2D  horizontal 
option  was  included  for  ease  of  use  and  allows  simplifying  assumptions,  but 
it  is  inappropriate  for  vertically  heterogeneous  datasets  of  columnar  points, 
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such  as  borehole  samples.  If  a  3D  interpolation  method  is  specified,  then 
full  3D  pilot  points  must  be  also  specified  via  PP  PT3.  The  3D  interpolation 
methods  (e.g.,  3D  ordinary  kriging)  require  the  use  of  PP  PT3. 

The  simplest  of  AdH’s  interpolation  types,  nearest-neighbor,  generates  a 
piecewise-constant  field  of  values  by  assigning  the  value  of  the  closest  pilot 
point.  The  minimum  and  maximum  number  of  pilot  points  (PP  RAD) 
must  be  equal  to  one  for  proper  specification  (enforced  to  avoid  any 
confusion).  Although  not  a  smooth  interpolant,  this  method  is  quick  and 
may  be  used  to  make  direct  substitution  of  values.  If  the  pilot  point  group 
contains  a  single  point,  and  the  search  radius  is  greater  than  the  domain 
extent’s  diagonal  measure  (i.e.,  the  search  radius  envelopes  the  entire 
domain),  the  value  is  automatically  stored  in  the  normal  material  property 
data  structure  instead  of  being  interpolated  and  stored  at  every  element. 

Inverse-distance  weighting,  one  of  the  simplest  linear  interpolation 
methods,  gives  more  significance  to  the  closer  pilot  points  and  less  to  the 
more  distant  by  using  the  distance  between  the  pilot  points  and  the  location 
of  interest.  AdH  modifies  Shepard’s  classical  formulation  by  setting  the 
arbitrary  weighting  exponent  to  2.0  and  including  the  partial  sample  R- 
sphere  mentioned  by  Franke  and  Nielson  (1980),  presented  here  as 
Equations  3  through  5. 

Modified  Shepard’s  method  IDW: 


1=1 


W 


i(Xloi)  = 


m 


iiXloi) 


mk  Xloi  >  = 


R-d(xloi,xk) 


R*d(xloi,xk) 


where: 

v  =  IDW  estimator 
xloi  =  location  of  interest 
Wj  =  IDW  weight 


(3) 

(4) 


(5) 
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v  =  observed  value 
Xj  =  pilot  point 
mb  nij  =  modified  weight 

R  =  R-sphere  radius  of  influence 
d(xioi,  xfc)=  distance  between  location  of  interest  and  pilot  point. 

The  R-sphere  radius  is  set  to  the  distance  of  the  farthest  pilot  point  in  the 
subset  to  remove  scaling  effects  from  the  computation.  The  downfall  of  the 
IDW  interpolant  is  that  it  produces  a  distribution  with  local  extrema  at  the 
pilot  points  and  values  trending  towards  the  mean  between  the  observa¬ 
tions,  which  may  not  necessarily  reproduce  the  understood  distribution. 

Kriging  is  a  set  of  geostatistical  methods  to  estimate  unknown  values 
based  on  variances  between  known  observations  (Equations  6-13). 
Ordinary  kriging  is  the  most  common  type  of  kriging  as  it  is  referred  to  as 
the  “best  linear  unbiased  estimator.”  “Best”  is  implied  here  “only  in  the 
least-squares  error  sense”  (Deutsch  and  Journel  1998)  for  minimizing  the 
error  variance,  erj.  “Ordinary  kriging  is  ‘linear’  because  its  estimates  are 
weighted  linear  combinations  of  the  available  data;  it  is  ‘unbiased’  since  it 
tries  to  have  mR,  the  mean  residual  or  error,  equal  to  o”  (Isaaks  and 
Srivastava  1989).  This  allows  ordinary  kriging  to  assume  an  unknown  and 
constant  mean  as  opposed  to  simple  kriging  where  the  mean  must  be 
specified,  which  is  difficult  with  data  collected  in  the  field. 

Kriging  general  equations: 


(6) 


1=1 


/in  a  n 

mR=-J2Ri  =  -^Vi-vi 


(7) 


°2R  =  z:J2(Ri-mRY 


n 


i—1 


(8) 


where: 

v  =  kriging  estimator 
xLoi  =  location  of  interest 
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Wj  =  kriging  weight 
v  =  observed  value 
Xj  =  pilot  point  location 
mR  =  mean  error 
Ri  =  error  {vt  -  v{) 
erj  =  error  variance. 

Ordinary  kriging  conditions: 


m„=0 

(9) 

!>.  =  1 

1=1 

(10) 

Ordinary  kriging  minimized  error  variance: 


2  2 
°R=°  ~ 


Yj»iCiM+V 


1=1 


—  o  —  w-D 


(id 


where: 

cr 2  =  variance 

Ci  ioi  =  covariance  of  pilot  point  and  location  of  interest 
=  Lagrange  parameter 
w  =  weights  vector 

D  =  covariance  matrix  (pilot  point  and  location  of  interest). 
Ordinary  kriging  system  of  equations: 


C-w  —  D 


(12) 


ci,i  -  c1>n  1 

m,' 

-  Cn,n  1 

= 

Cn;toi 

1  •••  1  0 

.  h  . 

1 

(13) 


where: 


C  =  covariance  matrix  (pilot  points) 
Ci  j  =  covariance  of  pilot  points. 


ERDC/CHL  TR-13-16 


21 


Ordinary  kriging  currently  is  AdH’s  most  computationally  expensive 
interpolation  scheme,  but  its  flexibility  creates  the  most  customizable 
distribution  with  sparse  data.  To  estimate  a  value,  the  covariances  between 
individual  pilot  point  locations  and  the  location  of  interest  are  computed. 
These  inform  the  linear  system  which  is  solved  for  the  weights.  Finally,  the 
influential  portions  of  the  observed  data,  based  on  the  weights,  are 
summed.  Additional  information  is  required  to  compute  the  covariances 
used  in  the  system  of  equations  (see  the  Kriging  Interpolation  Properties 
section).  For  detailed  discussion  on  ordinary  kriging,  see  Chapter  12  of 
Isaaks  and  Srivastava  (1989).  AdH’s  kriging  methods  are  based,  in  part,  on 
the  Geostatistical  Software  Library  (GSLIB)  (Deutsch  and  Journel  1998). 

When  3D  pilot  points  (PP  PT3)  are  provided,  it  is  fairly  easy  to  switch 
among  all  six  interpolation  methods,  compare  the  interpolants,  and 
estimate  the  model  sensitivity  to  the  interpolation  scheme. 

The  group  property  cards  (Table  5),  are  required  with  PP  PT2  or  PP  PT3 
(mutually  exclusive)  for  each  pilot  point  group. 

3. 1.2.3  Kriging  Interpolation  Properties 

These  cards  define  the  information  specific  to  the  kriging  interpolation 
scheme  and  are  necessary  only  if  the  PP  TYP  card  specifies  the  ordinary 
kriging  interpolation  type.  As  mentioned  in  the  general  interpolation 
discussion  above,  the  kriging  interpolant  estimates  values  influenced  by 
covariances  that  are  in  turn  derived  from  a  semivariogram  model.  A 
semivariogram,  or,  simply,  variogram,  is  a  description  of  spatial  variability 
as  a  function  of  distance,  y(h).  It  is  normally  presented  as  a  lD  scatter  plot 
of  observed  data  with  a  best-fit  curve.  For  purposes  in  AdH,  a  variogram 
must  be  viewed  as  a  3D  volume  or  at  least  a  2D  surface  description  to 
handle  anisotropy. 

The  key  parameters  of  a  variogram  are  as  follows: 

•  nugget  -  the  variance  value  at  h  =  0,  technically  should  be  equal  to 
zero  though  sampling  at  very  small  distances  may  cause  the  variance  to 
be  nonzero;  global  background  variance 

•  range  -  the  distance  at  which  the  change  in  variance  is  negligible 
(where  the  variogram  plateaus) 

•  sill  -  the  variance  value  at  the  range;  the  maximum  variance 
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•  bearing  angle  -  the  direction  of  the  major  axis  defined  by  the  three 
components,  measured  in  degrees: 

o  azimuth  -  rotation  around  the  z-axis  from  the  y-axis  (similar  to  yaw 
in  aeronautics)  where  clockwise  is  positive 
o  dip  -  rotation  around  the  x-axis  from  the  z-axis  (pitch)  where 
counterclockwise  is  positive 

o  plunge  -  rotation  around  the  y-axis  from  the  x-axis  (roll),  where 
clockwise  is  positive 

•  anisotropy  ratio  -  the  relationship  between  the  major  axis  and  two 
orthogonal  minor  axes: 

o  horizontal  -  the  non-rotated  y-axis  and  x-axis 
o  vertical  -  the  non-rotated  z-axis  and  y-axis. 


Table  5.  Group  Property  Cards 


PP  LIM 

INTERPOLATION  LIMITS 

The  PP  LIM  card  is  used  to  specify  the  interpolation  limits  for  a  pilot  point  group.  The 
interpolated  value  is  set  to  the  respective  limit  if  outside  the  given  valid  range.  If  interpolation 
cannot  be  performed  (e.g.,  search  criteria  are  not  met),  then  the  default  value  is  used. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

LIM 

Parameter 

3 

int 

*  INT_MAX 

Pilot  point  group  ID 

4 

int 

# 

Lower  interpolation  limit 

5 

int 

# 

Upper  interpolation  limit 

6 

int 

# 

Default  interpolation  value 

PPPT2 

2D  PILOT  POINTS 

The  PP  PT2  card  is  used  to  specify  a  group  of  2D  pilot  points,  each  with  a  location  and 
parameter  value.  The  number  of  lines  is  variable. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

PT2 

Parameter 

3 

int 

*  INT_MAX 

Pilot  point  group  ID 

4 

int 

>  1 

Number  of  points 

An  additional  line  with  the  following  fields  is  expected  for  each  point.  Blank  lines  or 
comment  lines,  demarcated  with  a  preceding  #  or  /,  may  be  interspersed,  but  the  total 
number  of  point  information  lines  must  be  provided  prior  to  any  other  card. 

Secondary  line  (for  each  point] 
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Field 

Type 

Value 

Description 

1 

string 

un 

Label  (use  single  or  double  quotes  for  multiple 
words] 

2 

real 

# 

X-coordinate 

3 

real 

# 

Y-coordinate 

4 

real 

# 

Value 

PPPT3 

3D  PILOT  POINTS 

The  PP  PT3  card  is  used  to  specify  a  group  of  3D  pilot  points,  each  with  a  location  and 
parameter  value.  The  number  of  lines  is  variable. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

PT3 

Parameter 

3 

int 

*  INT_MAX 

Pilot  point  group  ID 

4 

int 

>  1 

Number  of  points 

An  additional  line  with  the  following  fields  is  expected  for  each  point.  Blank  lines  or 
comment  lines,  demarcated  with  a  preceding  #  or  /,  may  be  interspersed,  but  the  total 
number  of  point  information  lines  must  be  provided  prior  to  any  other  card. 

Secondary  line  (for  each  point] 

Field 

Type 

Value 

Description 

1 

string 

im 

Label  (use  single  or  double  quotes  for  multiple 
words] 

2 

real 

# 

X-coordinate 

3 

real 

# 

Y-coordinate 

4 

real 

# 

Z-coordinate 

5 

real 

# 

Value 

PP  RAD 

SEARCH  RADIUS 

The  PP  RAD  card  is  used  to  specify  the  interpolation  search  parameters  for  a  pilot  point 
group.  These  parameters  define  the  subset  of  points  involved  in  the  interpolation  at  a  given 
domain  location.  If  more  points  than  the  maximum  (m)  are  found  within  the  search 
perimeter,  then  only  the  closest  m  points  are  involved. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

RAD 

Parameter 

3 

int 

t  INT_MAX 

Pilot  point  group  ID 

4 

real 

>0 

Search  Radius 

5 

int 

>  1 

Minimum  number  of  points 

6 

int 

>  minimum 

Maximum  number  of  points 

PPTYP 

INTERPOLATION  METHOD 

The  PP  TYP  card  is  used  to  specify  the  interpolation  method  for  a  pilot  point  group. 
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Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

TYP 

Parameter 

3 

int 

*  INT_MAX 

Pilot  point  group  ID 

4 

enum 

>0 

Interpolation  type: 

0  Nearest-neighbor 

I  Ordinary  kriging 

II  Modified  Shepard’s  method  inverse- 
distance  weighted  (IDW) 

5 

enum 

>0 

Interpolation  dimensionality: 

0  3D 

1  2D  horizontal 

An  omnidirectional  variogram,  defined  by  anisotropy  ratios  of  one,  depicts 
the  same  variance  in  any  direction  from  a  point.  In  this  case  the  specified 
bearing  angle  is  inconsequential.  When  the  ratios  are  non-uniform,  the 
bearing  angle  and  the  ratios  determine  a  transformation  matrix  that  scales 
and  rotates  the  variogram  to  compute  the  proper  y(/i).  It  may  be  easier  to 
understand  the  parameters  of  an  anisotropic  variogram  as  those  necessary 
to  transform  a  sphere  into  an  ellipsoid  with  semi-major  and  semi-minor 
axes  equal  to  the  range.  This  imaginary  ellipsoid  would  contain  the  volume 
within  which  the  variability  is  increasing. 

A  sample  variogram  is  usually  derived  during  the  analysis  of  observed  data 
showing  the  relationship  between  known  points.  For  kriging  purposes,  a 
variogram  is  necessary  to  describe  the  relationship  between  observed  data 
and  an  unobserved  location  in  the  determination  of  weights  where  the 
“variogram  distance  measures  the  average  degree  of  dissimilarity  between 
an  unsampled  value  z(it)  and  a  nearby  data  value.  For  example,  given  only 
two  data  values  z(it  +  h)  and  z(u  +  h ')  at  two  different  locations,  the  more 
dissimilar  sample  value  should  receive  lesser  weight  in  the  estimation  of 
z(it)”  (Deutsch  and  Journel  1998).  Used  in  this  fashion,  the  variogram  is 
referred  to  as  a  variogram  model. 

A  variogram  model  is  specified  using  a  best-fit  equation  to  express  the 
change  of  variance  between  the  nugget  and  sill.  For  complex  descriptions, 
multiple  expressions  can  be  summed,  referred  to  as  nested  variograms. 
The  available  variogram  types  are  presented  in  Equations  14  through  20. 
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Nugget-effect: 


Spherical: 


vih) 


\0,ifh  =  0 
\c,ifh>  0 


Y{h) 


h 

h 

3' 

c* 

1.5 - 0.5 

a 

a. 

c,ifh  >  0 

,ifh<  0 


Exponential: 


y(/z)  =  c* 


1-exp 


Gaussian: 


y(^)  — c 


c* 

1-exp 

W2  j 
2 

a 

{  )m 

Power: 


Y(h)  =  c*ha 


Hole  effect: 


y(/i)  =  c 


1  —  cos 

(h  ) 
—n 

a  ) 

(14) 


(15) 


(16) 


(17) 


(18) 


(19) 


Dampened  hole  effect: 


3  h) 

\h  11 

1-exp 

— 

-cos 

—  71 

d  J 

J 

y(/z)  =  c* 


(20) 
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where: 

y(h)  =  variogram  variance 
h  =  variogram  distance 
c  =  contribution 
a  =  range 

co  =  power  exponent,  0  <  a)  <  2 

d  =  distance  where  95%  of  the  hole  effect  is  dampened  out. 

The  contribution  (c)  of  the  variogram  type  is  the  same  concept  as  the  sill; 
however,  the  term  sill  is  reserved  here  for  only  the  variogram  model.  The  sill 
of  a  variogram  model  is  equal  to  the  sum  of  the  nugget  and  all  contributions. 
For  variogram  types  that  approach  their  contributions  asymptotically 
(exponential  and  Gaussian),  the  range  is  the  distance  at  which  y(a)  —  0.95c. 
The  power  variogram  type  does  not  technically  have  a  sill  (i.e.,  it  does  not 
plateau),  so  the  contribution  is  an  arbitrary  maximum  value.  An  anisotropic, 
nested  variogram  model  is  shown  in  Figure  9.  The  model  consists  of  a 
nugget  of  0.2;  an  exponential  type  variogram  with  a  contribution  of  0.3  and 
range  of  15;  and  a  Gaussian-type  variogram  with  a  contribution  of  0.5, 
major  axis  range  of  70,  horizontal  anisotropy  ratio  of  0.5,  and  azimuth 
bearing  angle  of  90  degrees.  The  major  axis  is  aligned  with  the  x-axis. 

To  provide  easy  and  flexible  specification  of  variogram  models,  parameters 
are  organized  on  multiple  cards.  An  individual  variogram’s  specific 
parameters  are  given  on  the  PP  VG2  or  PP  VG3  cards.  The  difference 
between  the  two  versions  is  the  user’s  approach  in  specifying  parameters; 
use  of  PP  VG2  assumes  that  vertical  transformation  and  rotations  are 
unnecessary,  and  AdH  defaults  them  appropriately.  The  PP  VG3  card 
requires  the  specification  of  all  three  bearing  angles  and  both  anisotropy 
ratios.  The  dimensionality  of  the  variogram  input  does  not  limit  the 
dimensionality  of  associated  kriging  interpolation.  Either  version  provides 
the  information  necessary  to  construct  a  variogram  to  describe  spatial 
variance  in  any  3D  direction.  The  first  characteristic  of  the  variogram 
model,  the  number  of  variograms  composing  the  model,  is  given  on  the  PP 
KRG  card.  Since  this  number  is  variable,  PP  VGI  lists  the  IDs  of  the 
variograms,  thereby  linking  the  pilot  point  group  to  PP  VG2  and  PP  VG3 
cards.  Second,  the  variogram  model  nugget  is  given  on  PP  KRG  which  is 
combined  with  the  third  characteristic,  the  individual  variogram’s 
contributions,  listed  on  the  PP  VGC  card,  to  produce  the  variogram  model 
sill.  The  contributions  must  be  listed  in  the  same  order  as  PP  VGI. 
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Figure  9.  Example  variogram  model  shown  as  a  surface  (A)  with  the  transformation  matrix  as 
an  ellipsoid  (B)  and  as  traditional  variogram  plots  along  the  principle  axes  (C  &  D). 


C)  Variogram  Model  Along  Major  Axis 


h 
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In  practice,  the  variance  magnitude  is  not  as  important  as  the  shape  of  the 
variogram.  Therefore,  the  sill  of  a  variogram  model  is  often  equal  to  one, 
y(a)  —  1,  for  the  sake  of  simplicity.  The  PP  VGW  card  is  provided  as  an 
alternative  to  (and  mutually  exclusive  with)  PP  VGC  for  those  inclined  to 
this  concept.  Instead  of  providing  contributions,  PP  VGW  lists  the 
weights  of  the  variograms  in  a  similar  fashion.  With  an  assumed  sill  and  a 
given  nugget,  the  weights  are  used  to  divide  the  remainder  and  compute 
contributions  for  the  user  as  defined  by  Equations  21  and  22. 

Variogram  contributions  from  weights: 

ci  =  A*[y(a)-y(0)]  (21) 

£ft=l  (22) 

1=1 


where: 

Cj  =  variogram  contribution 
Pi  =  variogram  weight 
y(a)  =  variogram  model  sill 
y(0)  =  variogram  model  nugget. 

The  weights  are  decimal  fractions  that  must  sum  to  1.0  (100%).  Due  to  the 
shape  of  the  expression,  PP  VGW  cannot  be  specified  when  including  a 
power-type  variogram.  If  necessary,  the  sill  assumption  can  be  overridden 
with  the  PP  VGS  card.  This  card  may  be  given  only  when  computing 
contributions  from  weights. 

As  stipulated  by  Isaaks  and  Srivastava  (1998),  “[i]f  the  major  features  of  the 
sample  variogram  can  be  captured  by  a  simple  model,  then  it  will  provide 
solutions  that  are  as  accurate  as  those  found  using  a  more  complex  model. 
The  principle  of  parsimony  is  a  good  guide  in  variogram  modeling.”  For 
detailed  discussion  on  modeling  variograms  see  Chapters  16  and  II. 3  of 
Isaaks  and  Srivastava  (1989)  and  Deutsch  and  Journel  (1998),  respectively. 

The  variogram  model  is  converted  to  a  covariance  function  for  the  purpose 
of  kriging  by  the  following  relationship,  Equation  23,  based  on  the 
assumption  that  the  random  variables  of  ordinary  kriging  all  have  the 
same  variance  and  mean. 
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Ordinary  kriging  variogram/covariance  relationship: 


cij  =  °2-  Y/j 


(23) 


where: 

Cij  =  covariance 
o2  =  kriging  variance 
Yij  =  variogram  model. 

Finally,  the  kriging  properties  also  include  an  option  on  the  PP  KRG  card 
to  utilize  a  log  transformation  in  the  interpolation.  This  may  provide  better 
results  when  interpolating  a  parameter  from  a  sample  with  a  range  of 
multiple  orders  of  magnitudes  (e.g.,  hydraulic  conductivity).  When 
selected,  the  general  equation,  Equation  6,  is  replaced  with  Equation  24. 

Kriging  equation  with  log  transformation: 


v[x, 


loi 


i)  =  10 


SU’i(Xtoi)loglo(u(xi)) 


(24) 


where: 

v  =  kriging  estimator 
xloi  =  location  of  interest 
Wj  =  kriging  weight 
v  =  observed  value 
X{  =  pilot  point  location. 

A  set  of  kriging  interpolation  property  cards  (Table  6)  is  required  for  each 
pilot  point  group  that  specifies  kriging  interpolation.  The  number  of 
variogram  cards,  however,  does  not  necessarily  correlate  the  number  of 
kriging  pilot  point  groups  as  variograms  may  be  nested  and/or  shared  by 
pilot  point  groups  (referenced  by  multiple  groups,  therefore  reusing  the 
same  spatial  variability).  PP  KRG  and  PP  VGI  are  required  with  one  of 
the  following  combinations:  PP  VGC;  PP  VGW;  or  PP  VGW  and  PP 
VGS.  Variogram  cards,  PP  VG2  and  PP  VG3,  are  mutually  exclusive  with 
at  least  one  required  to  perform  kriging. 
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Table  6.  Kriging  Interpolation  Property  Cards 


PP  KRG 

KRIGING  INFORMATION 

The  PP  KRG  card  is  used  to  specify  the  ordinary  kriging  interpolation  parameters  for  a  pilot 
point  group. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

KRG 

Parameter 

3 

int 

*  INT_MAX 

Pilot  point  group  ID 

4 

real 

>0 

Variogram  model  nugget 

5 

enum 

>0 

Transform: 

0  None 

1  Log 

6 

int 

>  1 

Number  of  variograms 

PPVGC 

VARIOGRAM  CONTRIBUTIONS 

The  PP  VGC  card  is  used  to  specify  the  contributions  for  variograms  used  in  kriging 
interpolation  for  a  pilot  point  group.  The  number  of  fields  is  variable.  The  number  of 
variogram  contributions  (n)  must  match  the  number  provided  by  the  PP  KRG  card.  The 
contributions  reference  the  variograms  listed  by  the  PP  VGI  card  (contributions  are  provided 
in  the  same  order  as  the  variogram  IDs). 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

VGC 

Parameter 

3 

int 

t  INT_MAX 

Pilot  point  group  ID 

4 

real 

>0 

1st  variogram  contribution 

3  +  n 

real 

>0 

nth  variogram  contribution 

PP  VGI 

VARIOGRAM  IDS 

The  PP  VGI  card  is  used  to  specify  which  variograms  are  used  for  the  kriging  interpolation 
for  a  pilot  point  group.  The  number  of  fields  is  variable.  The  number  of  variogram  IDs 
(n)must  match  the  number  provided  by  the  PP  KRG  card. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

VGI 

Parameter 

3 

int 

t  INT_MAX 

Pilot  point  group  ID 

4 

int 

*  INT_MAX 

1st  variogram  ID 

3  +  n 

int 

t  INT_MAX 

nth  variogram  ID 
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PPVGS 

VARIOGRAM  SILL 

The  optional  PP  VGS  is  used  to  specify  the  sill  of  the  variogram  model  used  in  kriging 
interpolation  for  a  pilot  point  group  when  contributions  are  computed  from  variogram 
weights. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

VGS 

Parameter 

3 

int 

t  INT_MAX 

Pilot  point  group  ID 

4 

real 

>0 

Variogram  model  sill 

PPVGW 

VARIOGRAM  WEIGHTS 

The  PP  VGW  card  is  used  to  specify  the  weights  for  variograms  used  in  kriging  interpolation  for 
a  pilot  point  group.  The  number  of  fields  is  variable.  The  number  of  variogram  weights  (n) 
must  match  the  number  provided  by  the  PP  KRG  card.  The  weights  reference  the  variograms 
listed  by  the  PP  VGI  card  (weights  are  provided  in  the  same  order  as  the  variogram  IDs). 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

VGW 

Parameter 

3 

int 

t  INT_MAX 

Pilot  point  group  ID 

4 

real 

>0 

1st  variogram  weight  (decimal  percent) 

3  +  n 

real 

>0 

nth  variogram  weight  (decimal  percent) 

PPVG2 

2D  VARIOGRAM  INFORMATION 

The  PP  VG2  card  is  used  to  specify  a  variogram  with  2D  parameters  for  use  with  kriging 
interpolation  (vertical  components  are  defaulted  appropriately).  The  last  fields  are 
conditional  on  a  preceding  field. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

VG2 

Parameter 

3 

int 

t  INT_MAX 

Pilot  point  group  ID 

4 

enum 

>0 

Variogram  type: 

0  Nugget-effect 

1  Spherical 

2  Exponential 

3  Gaussian 

4  Power 

5  Hole  effect 

6  Dampened  hole  effect 

5 

real 

>0 

Horizontal  anisotropy  ratio 

6 

real 

|#|  <360 

Azimuth  bearing  angle 

If  variogram  type  is  power  (field  4  is  equal  to  4]: 

7 

real 

0  <  #  <  2 

Variogram  power  exponent 
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Else: 

7 

real 

>0 

Major  axis  range 

If  variogram  type  is  dampened  hole  effect  (field  4  is  equal  to  6) 

8 

real 

>0 

Variogram  damping  lag  distance 

PP  VG3  3D  VARIOGRAM  INFORMATION 

The  PP  VG3  card  is  used  to  specify  a  variogram  with  3D  parameters  for  use  with  kriging 
interpolation.  The  last  fields  are  conditional  on  a  preceding  field. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

VG3 

Parameter 

3 

int 

t  INT_MAX 

Pilot  point  group  ID 

4 

enum 

>0 

Variogram  type: 

0  Nugget-effect 

1  Spherical 

2  Exponential 

3  Gaussian 

4  Power 

5  Hole  effect 

6  Dampened  hole  effect 

5 

real 

>0 

Horizontal  anisotropy  ratio 

6 

real 

>0 

Vertical  anisotropy  ratio 

7 

real 

|#|  <360 

Azimuth  bearing  angle 

8 

real 

|#|  <360 

Dip  bearing  angle 

9 

real 

|#|  <360 

Plunge  bearing  angle 

If  variogram  type  is  power  (field  4  is  equal  to  4]: 

10 

real 

0  <  #  <  2 

Variogram  power  exponent 

Else: 

10 

real 

>0 

Major  axis  range 

If  variogram  type  is  dampened  hole  effect  (field  4  is  equal  to  6) 

11 

real 

>0 

Variogram  damping  lag  distance 

3. 1.2.4  Miscellaneous 

The  following  PP  DBG  card,  listed  in  Table  7,  provides  supplemental 
information  regarding  pilot  point  specification.  It  is  unlikely  to  be  used  in 
general  practice  as  its  intent  is  to  summarize  input  for  testing  purposes. 
For  the  sake  of  completeness,  the  card  is  included  in  this  documentation. 
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Table  7.  Miscellaneous  Cards 


PP  DBG 

DEBUG  INFORMATION 

The  presence  of  the  optional  PP  DBG  card  causes  AdH  to  print  additional  pilot  point 
information  with  the  regular  screen  output  for  review  of  input  and  to  assist  debugging. 

Field 

Type 

Value 

Description 

1 

string 

PP 

Card  type 

2 

string 

DBG 

Parameter 

3.2  Output 

AdH  will  output  the  interpolated  pilot  point  parameter  values  as  datasets 
for  inspection  and  verification  in  the  eXentsible  Data  Model  and  Format 
(XDMF1,  maintained  by  Kitware,  Inc.)  file  format  that  is  viewable  in  the 
ParaView2 3  visualization  software  (developed  by  Kitware,  Inc.;  Henderson 
2007).  The  XDMF  file  format  is  a  combination  of  XML  light  data  and 
HDF5  heavy  data  so  the  raw  data  arrays  can  also  be  inspected  with  the 
HDFView3  tool  (developed  by  The  HDF  Group).  Since  mesh  element  data 
sets  are  not  supported  by  the  compatible  ASCII  file  format,  AdH  cannot 
use  the  traditional  output  to  provide  pilot  point  data  to  the  Department  of 
Defense  Groundwater  Modeling  System  (GMS,  developed  by  Aquaveo, 
LLC). 

The  resulting  solution  datasets  of  an  AdH  simulation  (e.g.,  total  head)  are 
written  as  normal,  incorporating  the  effects  of  the  pilot  point  specification. 


1  http://www.xdmf.org/ 

2  http://www.paraview.org/ 

3  http://www.hdfgroup.org/HDF5/ 
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4  Running  AdH  with  Pilot  Points 

When  all  required  input  is  ready,  an  AdH  simulation  utilizing  pilot  point 
specification  is  executed  in  the  normal  fashion  with  a  pilot  point-enabled 
version  of  the  code.  To  verify  the  version  of  AdH  executable,  call  the 
executable  with  the  argument  -v  as  shown  with  the  resulting  build 
information  in  Figure  to.  The  PILOT_POINTS  keyword  will  be  listed  if 
enabled.  Likewise,  the  inclusion  of  pilot  point  interpolated  parameter 
fields  in  the  output  is  dependent  on  the  XDMF  keyword. 

Figure  10.  Verifying  AdH  executable  is  pilot  point  enabled 
in  a  UNIX  shell  (Bash). 

$  . /adh  -v 

AdH  Build  Information 
SVN  revision  #  12324 

Build  Date/Time:  2013.06.26  /  14:24:20 
Built  with  GW  physics  enabled 
Built  with  XDMF  output  file  format 
Built  with  MPI  enabled 
Built  with  PILOT_POINTS  enabled 


Although  not  necessary  for  regular  groundwater  problems,  the  pre-AdH 
executable  may  be  run  first  to  verify  interpolated  parameters.  The 
difference  between  the  pre-AdH  and  AdH  executables  is  that  the  former 
does  not  proceed  into  the  computational  loops.  Both  executables  read  and 
verify  the  input,  initialize  the  problem,  perform  pilot  point  interpolation, 
and  write  the  interpolated  parameters  to  the  output.  Pre-AdH  may  be 
verified  and  executed  in  a  similar  fashion  to  the  full  executable. 

To  run  AdH,  call  the  executable  with  the  simulation  base  name  as  an 
argument  (Figure  u).  The  executable  build  information  is  written  to  the 
screen,  followed  by  runtime  information,  simulation  input  information, 
initialization,  and,  finally,  the  computational  proceedings  as  shown  in 
Figures  n  and  12.  The  former  shows  the  head  of  the  AdH  output,  while  the 
latter  shows  the  continuation  of  the  output,  though  truncated,  including  the 
beginning  of  the  computational  loops.  Check  the  runtime  information 
section  to  ensure  the  proper  simulation  input  files  were  provided  to  AdH. 
Either  default  or  specified  filenames  will  be  listed,  depending  on  whether  an 
AdH  super  file  is  found.  These  are  the  files  AdH  will  read  to  retrieve 
necessary  input.  Pilot  point  files  will  only  be  read  if  AdH  is  directed  to 
include  pilot  point  specification. 
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Figure  11.  Running  AdH  executable  in  a  UNIX  shell  (Bash). 


$  . /adh  Ex_sim 


AdH  Build  Information 


SVN  revision  #  12324 

Build  Date/Time:  2013.06.26  /  14:24:20 
Built  with  GW  physics  enabled 
Built  with  XDMF  output  file  format 
Built  with  MPI  enabled 
Built  with  PILOT_POINTS  enabled 


Runtime  Information 


AdH  execution  Date/Time:  2013.05.23  /  14:30:37 

Launching  AdH  with  project  name:  Ex_sim 

Launching  AdH  with  run  name:  Ex_sim 

Found  Super  file  named:  Ex_sim.sup 

Default  Geometry  file  name:  Ex_sim.3dm 

Default  Boundary  condition  file  name:  Ex_sim.bc 

Default  Groundwater  hots tart  file  name:  Ex_sim.hot 

Specified  Pilot  Point  file  name:  K-pts . txt  (used  if  OP  PP  is  specified) 
Specified  Pilot  Point  file  name:  K-Kriged.txt  (used  if  OP  PP  is  specified) 
AdH  was  launched  with  1  processor 


Geometry  Information 


Number  of  2D  elements:  0 
Number  of  3D  elements:  54000 
Number  of  nodes:  10461 


At  runtime,  AdH  will  read  the  pilot  point  files  only  after  reading  the 
standard  input  and  finding  the  OP  PP  card.  The  PP  cards  are  read  and 
validated  individually  to  confirm  that  card-level  input  criteria  have  been 
met.  Then,  each  specified  pilot  point  group  is  validated  to  ensure  all 
specifications  are  complete  and  coherent.  Any  issues  found  during  the 
validation  routine  will  be  listed  in  the  Pilot  Point  Information  section  of 
the  screen  output  (Figure  12)  and  cause  AdH  to  terminate  prematurely. 
Warnings  regarding  pilot  point  input,  which  do  not  require  resolution,  are 
also  listed  there.  AdH  attempts  to  validate  all  pilot  point  input  in  a  single 
pass  thus  reducing  the  cycle  of  fixing  one  error  only  to  find  yet  another 
during  the  subsequent  model  call.  When  all  pilot  point  input  is  acceptable, 
a  summary  of  the  pilot  point  specification  is  provided  in  the  Pilot  Point 
Information  section. 

Next,  AdH  will  link  the  pilot  point  groups  to  model  parameters  while 
assessing  compatibility  and  uniqueness  during  the  initialization  routine. 
Compatibility  includes  pilot  point  group-to-pilot  point  group  relationships 
(e.g.,  hydraulic  conductivity  components)  and  pilot  point  group-to-model 
relationships  such  as  the  existence  of  referenced  materials  and  direct 


ERDC/CHL  TR-13-16 


36 


Figure  12.  Example  AdH  pilot  point  information  screen  output  truncated  with  ellipses. 


Pilot  Point  Information 


Reading  file:  K-pts . txt 
Reading  file:  K-Kriged.txt 
Pilot  point  group  200: 

The  kriging  variogram  model  sill  is  assumed  to  be  1.0. 

Number  of  pilot  point  groups:  3 
Specified  interpolated  parameters: 

Type:  Material  property 

Material  ID:  1  Parameter:  Hydraulic  conductivity  scaling  factor 

Material  ID:  2  Parameter:  Hydraulic  conductivity  scaling  factor 

Material  ID:  3  Parameter:  Hydraulic  conductivity  scaling  factor 

Specified  interpolation  methods: 

2D  Horizontal  Ordinary  Kriging 
Total  number  of  pilot  points:  15 
Number  of  variograms :  3 


Pilot  Point  Initialization 


Initializing  and  linking  pilot  point  data. 

Computing  the  point  pair  covariances  for  pilot  point  group  100. 

Computing  the  point  pair  covariances  for  pilot  point  group  200. 

Computing  the  point  pair  covariances  for  pilot  point  group  300. 

Interpolating  the  3D  elements'  pilot  point  parameters. 

Completed  pilot  point  data  initialization. 

Printing  solution  at  time:  0.00000e+00  0.00000e+00  Percent  Done 


The  Master  Time  Loop 


************************************** 

*  Time  Interval  (0.000000,  1.000000)  * 
************************************** 


substitution  with  nearest-neighbor  interpolation  type.  Once  again,  if  any 
issues  are  found,  they  will  be  reported  to  the  screen,  this  time  in  the  Pilot 
Point  Initialization  section  (Figure  12),  and  AdH  will  exit.  Each  pilot  point 
group  is  interpolated  to  its  respective  domain  entity  at  this  time.  Computed 
initial  condition  values  are  assigned  directly  to  those  mesh  nodes  without 
existing  boundary  condition  assignments  (which  were  already  applied  from 
standard  input)  for  the  initial  solve.  Interpolated  material  property  values 
are  stored  in  an  auxiliary  array  structure  that  is  linked  to  by  the  material 
data  structure.  If  kriging  is  performed  during  the  initial  interpolation,  the 
covariances  between  pilot  point  pairs  are  computed  and  also  stored  for 
future  use. 

Prior  to  the  computational  loop,  AdH  writes  the  model  state  of  initial  and 
boundary-condition  forcings  as  the  first  solution  output.  The  initial 
condition  values  based  on  pilot  point  specification  are  integrated  into  this 
first  time-step  of  the  data  sets.  All  material  property  pilot  point-based 
values  are  also  written  at  this  time  as  element  field  data  with  NaN  (Not  a 
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Number)  reported  at  elements  utilizing  the  standard  zonal  material 
property  value. 

As  AdH  sets  up  and  solves  the  groundwater  or  heat  transport  problem,  it 
normally  retrieves  material  property  parameters  from  the  standard  data 
structure.  However,  when  pilot  point  specification  is  prescribed,  AdH  first 
checks  whether  a  particular  material  property  is  linked  to  a  spatially- 
varied  data  set  and  then  retrieves  the  value  from  the  auxiliary  pilot  point 
array  structure  or  the  standard  data  structure  as  necessary.  As  expected, 
this  retrieval  process  taxes  the  simulation  speed  slightly,  but  the  overall 
pilot  point  process  was  designed  to  balance  speed  with  memory  usage. 
Pilot  point  specification  is  integrated  into  AdH’s  domain  decomposition 
for  parallel  processing  and  domain  adaption  routines.  Every  processor 
holds  the  pilot  point  specification  input  information,  but  the  spatially- 
varying  data  held  by  the  auxiliary  array  structure  is  limited  to  each 
processor’s  subdomain.  When  the  mesh  is  refined  or  unrefined,  the 
auxiliary  array  structure  is  updated  to  match  the  number  of  elements  and 
new  pilot  point-based  values  are  interpolated  for  each  adapted  element 
since  their  centroids  changed.  The  new  pilot  point  data  is  also  written 
when  the  new  mesh  topology  is  outputted. 

Given  the  values  assigned  through  the  pilot  point  process,  AdH  will 
operate  and  finish  normally.  Pilot  point  specification  only  substitutes 
values;  the  use  of  these  values  is  not  altered. 
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Appendix:  Pilot  Point  Input  File  Examples 

Example  A 


Simulation  Ex_sim_A  shows  pilot  point  specification  in  default  location 
and  consists  of  the  following  files: 


•  Ex__sim_A.3dm  -  mesh  geometry  file 

•  Ex__sim_A.bc  -  boundary  condition  file  (shown  below) 

•  Ex_sim_A.hot  -  hotstart,  or  initial  condition,  file 

•  Ex_sim_A.pp  -  pilot  point  specification  file  (shown  below) 

Ex_sim_A.bc  (truncated  with  ellipses) 


OP  GW 

OP  PP  #  Enable  pilot  point  specification 
OP  TRN  0 

MP  G  73231257600.0  #  Gravity,  m/dayA2  (9.81  m/sA2) 

MP  RHO  1060.0  #  Density,  kg/mA3  (1.06e6  g/mA3) 

MP  MU  86.5728  #  Viscosity,  g/(m*s)  (1.002  g/(m*s) 

#  Material  1  -  Top  Layer 

MP  K  1  3.0e-l  3.0e-l  3 . Oe-2  0.0  0.0  0.0  #  Hydraulic  conductivity,  m/day 

MP  POR  10.4#  Porosity 

MP  SS  1  1.0e-5  #  Specific  storage,  1/m 


Ex_sim_A.pp 


#  Material  1  Hydraulic  conductivity  scaling 

PP  MP  KS  100  1  #  parameter,  pilot  point  group  id,  material  id 
PP  TYP  100  11  1  #  ppg  id,  interp  type  (IDW) ,  interp  dim  (2D  horiz) 

PP  RAD  100  120.0  25#  ppg  id,  search  radius,  min  pts,  max  pts 

PP  LIM  100  0.1  10.0  -999  #  ppg  id,  low  bound,  high  bound,  default  value 

PP  PT2  100  5  #  ppg  id,  num  pts 

SW  0.0  0.0  2.15  #  pt  label,  X,  Y,  value 

SE  90.0  0.0  3.37 

Well  60.0  30.0  1.00 

NE  90.0  60.0  0.92 

NW  0.0  60.0  0.22 


Example  B 

Simulation  Ex_sim_B  shows  pilot  point  specification  in  the  ADH  input 
file,  referenced  by  the  super  file,  and  consists  of  the  following  files: 


•  Ex_sim_B.3dm  -  mesh  geometry  file 

•  Ex_sim_B.bc  -  boundary  condition  file  (shown  below) 

•  Ex_sim_B.hot  -  hotstart,  or  initial  condition,  file 
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•  Ex_sim_B.sup  -  simulation  super  file  (shown  below) 


Ex_sim_B.bc  (truncated  with  ellipses) 


OP  GW 

OP  PP  #  Enable  pilot  point  specification 
OP  TRN  0 

MP  G  73231257600.0  #  Gravity,  m/dayA2  (9.81  m/sA2) 

MP  RHO  1060.0  #  Density,  kg/mA3  (1.06e6  g/mA3) 

MP  MU  86.5728  #  Viscosity,  g/(m*s)  (1.002  g/(m*s) 

#  Material  1  -  Top  Layer 

MP  K  1  3.0e-l  3.0e-l  3 . Oe-2  0.0  0.0  0.0  #  Hydraulic  conductivity,  m/day 

MP  POR  1  0.4  #  Porosity 

MP  SS  1  1.0e-5  #  Specific  storage,  1/m 

#  Material  1  Hydraulic  conductivity  scaling 

PP  MP  KS  100  1  #  parameter,  pilot  point  group  id,  material  id 
PP  TYP  100  11  1  #  ppg  id,  interp  type  (IDW) ,  interp  dim  (2D  horiz) 

PP  RAD  100  120.0  25#  ppg  id,  search  radius,  min  pts,  max  pts 

PP  LIM  100  0.1  10.0  -999  #  ppg  id,  low  bound,  high  bound,  default  value 

PP  PT2  100  5  #  ppg  id,  num  pts 

SW  0.0  0.0  2.15  #  pt  label,  X,  Y,  value 

SE  90.0  0.0  3.37 

Well  60.0  30.0  1.00 

NE  90.0  60.0  0.92 

NW  0.0  60.0  0.22 


Ex_sim_B.sup 


PP  Ex_s im_B . be 


Example  C 

Simulation  Ex_sim_C  shows  pilot  point  specification  in  multiple  files, 
referenced  by  the  super  file,  and  consists  of  the  following  files: 

•  Ex_sim_C.3dm  -  mesh  geometry  file 

•  Ex_sim_C.bc  -  boundary  condition  file  (shown  below) 

•  Ex_sim_C.hot  -  hotstart,  or  initial  condition,  file 

•  Ex_sim_C.sup  -  simulation  super  file  (shown  below) 

•  Mati_K_IDW.txt  -  pilot  point  specification  file  (shown  below) 

•  Mati_K_pts.txt  -  pilot  point  specification  file  (shown  below) 


ERDC/CHL  TR-13-16 


41 


Ex_sim_C.bc  (truncated  with  ellipses) 


OP  GW 

OP  PP  #  Enable  pilot  point  specification 
OP  TRN  0 

MP  G  73231257600.0  #  Gravity,  m/dayA2  (9.81  m/sA2) 

MP  RHO  1060.0  #  Density,  kg/mA3  (1.06e6  g/mA3) 

MP  MU  86.5728  #  Viscosity,  g/(m*s)  (1.002  g/(m*s) 

#  Material  1  -  Top  Layer 

MP  K  1  3.0e-l  3.0e-l  3 . Oe-2  0.0  0.0  0.0  #  Hydraulic  conductivity,  m/day 

MP  POR  10.4#  Porosity 

MP  SS  1  1.0e-5  #  Specific  storage,  1/m 


Ex_sim_  C.sup 


PP  Matl_K_pts . txt 
PP  Matl_K_IDW.txt 


Mail  K  IDW.txt 


#  Material  1  Hydraulic  conductivity  scaling 

PP  MP  KS  100  1  #  parameter,  pilot  point  group  id,  material  id 
PP  TYP  100  11  1  #  ppg  id,  interp  type  (IDW) ,  interp  dim  (2D  horiz) 

PP  RAD  100  120.0  25#  ppg  id,  search  radius,  min  pts,  max  pts 

PP  LIM  100  0.1  10.0  -999  #  ppg  id,  low  bound,  high  bound,  default  value 


Mati_K_pts.  txt 


#  Material  1  Hydraulic  conductivity  scaling 

PP  PT2  100  5  #  ppg  id,  num  pts 

SW  0.0  0.0  2.15  #  pt  label,  X,  Y,  value 

SE  90.0  0.0  3.37 

Well  60.0  30.0  1.00 

NE  90.0  60.0  0.92 

NW  0.0  60.0  0.22 


Example  D 

Simulation  Ex_sim_D  shows  ordinary  kriging  interpolation  method  pilot 
point  specification  in  the  default  location  and  consists  of  the  following 
files: 

•  Ex_sim_D.3dm  -  mesh  geometry  file 

•  Ex_sim_D.bc  -  boundary  condition  file  (shown  below) 

•  Ex_sim_D.hot  -  hotstart,  or  initial  condition,  file 

•  Ex_sim_D.pp  -  pilot  point  specification  file  (shown  below) 
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Ex_sim_D.bc  (truncated  with  ellipses) 


OP  GW 

OP  PP  #  Enable  pilot  point  specification 
OP  TRN  0 

MP  G  73231257600.0  #  Gravity,  m/dayA2  (9.81  m/sA2) 

MP  RHO  1060.0  #  Density,  kg/mA3  (1.06e6  g/mA3) 

MP  MU  86.5728  #  Viscosity,  g/(m*s)  (1.002  g/(m*s) 

#  Material  1  -  Top  Layer 

MP  K  1  3.0e-l  3.0e-l  3 . Oe-2  0.0  0.0  0.0  #  Hydraulic  conductivity,  m/day 

MP  POR  10.4#  Porosity 

MP  SS  1  1.0e-5  #  Specific  storage,  1/m 


Ex_sim_D.pp 


#  Variograms 

#  variogram  id,  variogram  type  (gsn) ,  horiz  anisotropy,  bearing  angle,  range 
PP  VG2  -100  3  1.0  0.0  15.0 

PP  VG2  -200  2  0.5  90.0  70.0 

#  Material  1  Hydraulic  conductivity  scaling 

PP  MP  KS  100  1  #  parameter,  pilot  point  group  id,  material  id 

PP  TYP  100  11#  ppg  id,  interp  type  (kriging) ,  interp  dim  (2D  horiz) 

PP  KRG  100  0.2  02#  ppg  id,  nugget,  transform  (log),  num  variograms 
PP  VGI  100  -100  -200  #  ppg  id,  1st  variogram  id,  . . . 

PP  VGC  100  0.3  0.5  #  ppg  id,  1st  variogram  contrib,  ... 

PP  RAD  100  120.0  25#  ppg  id,  search  radius,  min  pts,  max  pts 

PP  LIM  100  0.1  10.0  -999  #  ppg  id,  low  bound,  high  bound,  default  value 

PP  PT2  100  5  #  ppg  id,  num  pts 

SW  0.0  0.0  2.15  #  pt  label,  X,  Y,  value 

SE  90.0  0.0  3.37 

Well  60.0  30.0  1.00 

NE  90.0  60.0  0.92 

NW  0.0  60.0  0.22 
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